###### Title ##########################
## Insecurity and Support for Female Leadership in Conflict States: 
## Evidence from Afghanistan by Jasmine Bhatia and
## Steve L. Monroe.
## Replication Command File: Conjoint Experiment  
## August 2023
##

# This file contains analysis and tables of the conjoint survey experiment
# presented in both the main article and supplementary information


### Table of Contents  ##################################
## 1. Load packages, data and clean data for MM analysis
## 2. Format data for Conjoint Analysis
## 3. H1 Respondents with the Insecurity Prime
##     will be less supportive of female leaders
##     A. H1 Table and Plot in Appendix
##     B. Power Analysis
## 4. H2 Insecurity Prime Across Genders
##     A. H2: Table across gender and treatment groups
##     B. H2a: Men who receive the prime are less supportive of female leaders 
##        a. Table and Plot
##        b. Power Analysis  
##     C. H2b: Women who receive the prime are less supportive of female leader
##       a. Table and Plot in the Appendix
##       b. Power Analysis
## 5. H2 Robustness Checks Subsets
##     a. Unemployed (F test, Table, Male and Female Respondents) 
##     b, Less Educated (F test, Table, Male and Female Respondents) 
##     c. Not living in Kundoz (F test, Table, Male and Feamale Respondents) 




## Note: All code run using R 4.0.3 GUI 1.73 Catalina build (7892)

######## 1. Load Necessary Packages and Data ##############################

library(cjoint)
library(dplyr)
library(tidyr)
library(xtable)
library(tidyverse)
library(stats)
library(survey)
library(Rcpp)
library(estimatr)
library(knitr)
library(DT)
library(cli)
library(devtools)
devtools::install_github("m-freitag/cjpowR")
library(cjpowR)
library(stargazer)
library(MASS)
library(mlogit)
library("nnet")
library("scales")
library(here)

# set working directory and load dataset

getwd()


data <- read.csv("AfghanDataInscCtrlOnly_Clean.csv")


#### 2. Format Data for Conjoint Analysis #####################

# Create separate lists of attributes for each 
# leadership profile

profile_1 <- c("Rd_1_A_Gender", "Rd_1_A_Education", "Rd_1_A_Age", 
               "Rd_1_A_Ethnicity", "Rd_1_A_Professional_Experience",
               "Rd_1_A_Place_of_Birth")

profile_2 <- c("Rd_1_B_Gender", "Rd_1_B_Education", "Rd_1_B_Age", 
               "Rd_1_B_Ethnicity", "Rd_1_B_Professional_Experience",
               "Rd_1_B_Place_of_Birth")
  
profile_3 <- c("Rd_2_A_Gender", "Rd_2_A_Education", "Rd_2_A_Age", 
               "Rd_2_A_Ethnicity", "Rd_2_A_Professional_Experience",
               "Rd_2_A_Place_of_Birth")

profile_4 <- c("Rd_2_B_Gender", "Rd_2_B_Education", "Rd_2_B_Age", 
               "Rd_2_B_Ethnicity", "Rd_2_B_Professional_Experience",
               "Rd_2_B_Place_of_Birth")

profile_5 <- c("Rd_3_A_Gender", "Rd_3_A_Education", "Rd_3_A_Age", 
               "Rd_3_A_Ethnicity", "Rd_3_A_Professional_Experience",
               "Rd_3_A_Place_of_Birth")

profile_6 <- c("Rd_3_B_Gender", "Rd_3_B_Education", "Rd_3_B_Age", 
               "Rd_3_B_Ethnicity", "Rd_3_B_Professional_Experience",
               "Rd_3_B_Place_of_Birth")

 
# Create separate list of responses for each pair of profiles

responses_1 <- c("RD1_Con_Table_Choice", "RD_1_A_Rank", "RD_1_B_Rank")
responses_2 <- c("RD2_Con_Table_Choice", "RD_2_A_Rank", "RD_2_B_Rank")  
responses_3 <- c("RD3_Con_Table_Choice", "RD_3_A_Rank", "RD_3_B_Rank")    


# Subset relevant treatments, responses, and all covariates 
# for each profile
covars <- c("ControlGrp", "TreatInsc", "Res_Gender", "Literacy", "HHhead",
            "Employ_Unem", "District", "HH_Head_Edu", "Res_Education",
            "Province", "Income", "Voted_Pres", "Res_Ethnicity", "Behav_Peace",
            "IntForces", "TrustProvGov", "TrustCentGov") 
  
data2_1 <- data[, c("ID", "Rd_1_A_Gender", "Rd_1_A_Education", "Rd_1_A_Age", 
                    "Rd_1_A_Ethnicity", "Rd_1_A_Professional_Experience",
                    "Rd_1_A_Place_of_Birth",
                    "RD1_Con_Table_Choice", "RD_1_A_Rank", "RD_1_B_Rank",
                    "ControlGrp", "TreatInsc", "Res_Gender", "Literacy", "HHhead",
                     "Employ_Unem", "District", "HH_Head_Edu", "Res_Education",
                    "Province", "Income", "Voted_Pres", "Res_Ethnicity",
                    "Behav_Peace", "IntForces", "TrustProvGov", "TrustCentGov")]
                    

data2_1$profile <- rep(1, 1299)                    
                    

data2_1$Pair <- "A" #name of pair. Will use this to create a profile pair id

data2_2 <- data[, c("ID", "Rd_1_B_Gender", "Rd_1_B_Education", "Rd_1_B_Age", 
                    "Rd_1_B_Ethnicity", "Rd_1_B_Professional_Experience",
                    "Rd_1_B_Place_of_Birth",
                    "RD1_Con_Table_Choice", "RD_1_A_Rank", "RD_1_B_Rank",
                    "ControlGrp", "TreatInsc", "Res_Gender", "Literacy", "HHhead",
                    "Employ_Unem", "District", "HH_Head_Edu", "Res_Education",
                    "Province", "Income", "Voted_Pres", "Res_Ethnicity",
                    "Behav_Peace", "IntForces", "TrustProvGov", "TrustCentGov")]


data2_2$profile <- rep(2, 1299)   

data2_2$Pair <- "A" #name of pair. Will use this to create a profile pair id

data2_3 <- data[, c("ID", "Rd_2_A_Gender", "Rd_2_A_Education", "Rd_2_A_Age", 
                    "Rd_2_A_Ethnicity", "Rd_2_A_Professional_Experience",
                    "Rd_2_A_Place_of_Birth",
                    "RD2_Con_Table_Choice", "RD_2_A_Rank", "RD_2_B_Rank",
                    "ControlGrp", "TreatInsc", "Res_Gender", "Literacy", "HHhead",
                    "Employ_Unem", "District", "HH_Head_Edu", "Res_Education",
                    "Province", "Income", "Voted_Pres", "Res_Ethnicity",
                    "Behav_Peace", "IntForces", "TrustProvGov", "TrustCentGov")]


data2_3$profile <- rep(3, 1299)                    

data2_3$Pair <- "B" #name of pair. Will use this to create a profile pair id

data2_4 <- data[, c("ID", "Rd_2_B_Gender", "Rd_2_B_Education", "Rd_2_B_Age", 
                    "Rd_2_B_Ethnicity", "Rd_2_B_Professional_Experience",
                    "Rd_2_B_Place_of_Birth",
                    "RD2_Con_Table_Choice", "RD_2_A_Rank", "RD_2_B_Rank",
                    "ControlGrp", "TreatInsc", "Res_Gender", "Literacy", "HHhead",
                    "Employ_Unem", "District", "HH_Head_Edu", "Res_Education",
                    "Province", "Income", "Voted_Pres", "Res_Ethnicity",
                    "Behav_Peace", "IntForces", "TrustProvGov", "TrustCentGov")]

data2_4$profile <- rep(4, 1299)  

data2_4$Pair <- "B" #name of pair. Will use this to create a profile pair id

data2_5 <- data[, c("ID", "Rd_3_A_Gender", "Rd_3_A_Education", "Rd_3_A_Age", 
                    "Rd_3_A_Ethnicity", "Rd_3_A_Professional_Experience",
                    "Rd_3_A_Place_of_Birth",
                    "RD3_Con_Table_Choice", "RD_3_A_Rank", "RD_3_B_Rank",
                    "ControlGrp", "TreatInsc", "Res_Gender", "Literacy", "HHhead",
                    "Employ_Unem", "District", "HH_Head_Edu", "Res_Education",
                    "Province", "Income", "Voted_Pres", "Res_Ethnicity",
                    "Behav_Peace", "IntForces", "TrustProvGov", "TrustCentGov")]


data2_5$profile <- rep(5, 1299)                    

data2_5$Pair <- "C" #name of pair. Will use this to create a profile pair id

data2_6 <- data[, c("ID", "Rd_3_B_Gender", "Rd_3_B_Education", "Rd_3_B_Age", 
                    "Rd_3_B_Ethnicity", "Rd_3_B_Professional_Experience",
                    "Rd_3_B_Place_of_Birth",
                    "RD3_Con_Table_Choice", "RD_3_A_Rank", "RD_3_B_Rank",
                    "ControlGrp", "TreatInsc", "Res_Gender", "Literacy", "HHhead",
                    "Employ_Unem", "District", "HH_Head_Edu", "Res_Education",
                    "Province", "Income", "Voted_Pres", "Res_Ethnicity",
                    "Behav_Peace", "IntForces", "TrustProvGov", "TrustCentGov")]


data2_6$profile <- rep(6, 1299)  

data2_6$Pair <- "C" #name of pair. Will use this to create a profile pair id

# Rename columns and rbind datasets back together
neutral_names <- c("ID", "Gender", "Education", "Age", 
                   "Ethnicity", "Career", "POB",
                   str_replace(responses_1, "_1", ""), 
                   covars, "profile", "Pair")

names(data2_1) <- neutral_names
names(data2_2) <- neutral_names
names(data2_3) <- neutral_names
names(data2_4) <- neutral_names
names(data2_5) <- neutral_names
names(data2_6) <- neutral_names


data2 <- data.frame(rbind(data2_1, data2_2, 
                          data2_3, data2_4, 
                          data2_5, data2_6 
                          ))

# Add a unique row identifier
data2$num <- 1:7794

## Create Pair Variable

data2$PairID <- paste(data2$ID, data2$Pair, sep="")

## Put variables in correct format for cjoint package functions

# Rename data for easier typing: 
data <- data2 


# Consolidate rating variables
data$rating <- case_when(
  data$profile %in% c(1,3,5) ~ data$RD_A_Rank,
  data$profile %in% c(2,4,6) ~ data$RD_B_Rank
)


# Note: NAs must be coded as zero for cjoint functions to work

data$choice_dummy <- ifelse(
        data$profile %in% c(1,3,5) & data$RD1_Con_Table_Choice == 1 |
        data$profile %in% c(2,4,6) & data$RD1_Con_Table_Choice == 2, 1, 0)


## Make Control Variables

# Military Variable Attribute

data$Military <- ifelse(data$Career == "Military", "Military", "Civilian")
data$Military <- as.factor(data$Military)


##  Gender Variable Attribute (make it a factor)

data$Gender <- as.factor(data$Gender)


# Unemployed Variable (Respondent)

# unemployed = 1
data$Employ_Unem[is.na(data$Employ_Unem)] <- 0


# Education level variable (Respondent)

data$edu_numeric <- ifelse(data$Res_Education == 1, 0,
                       ifelse(data$Res_Education == 8, 0,
                       ifelse(data$Res_Education == 2, 1,
                        ifelse(data$Res_Education == 3, 2,      
                        ifelse(data$Res_Education == 4, 3,     
                        ifelse(data$Res_Education == 5, 4,
                         ifelse(data$Res_Education == 6, 5, 
                         NA)))))))

# Make a new peace variable

data$Behav_Peace_new <- case_when(
  data$Behav_Peace == 1 ~ 6,
  data$Behav_Peace  == 2 ~ 5, 
  data$Behav_Peace == 3 ~ 4, 
  data$Behav_Peace == 4 ~ 3,
  data$Behav_Peace == 5 ~ 2,
  data$Behav_Peace== 6 ~ 1,
  TRUE ~ NA_real_
)


# Make a new international forces support variable

data$IntForces_new <- case_when(
  data$IntForces == 1 ~ 5,
  data$IntForces  == 2 ~ 4, 
  data$IntForces == 3 ~ 3, 
  data$IntForces == 4 ~ 2,
  data$IntForces == 5 ~ 1,
  TRUE ~ NA_real_
)


# Make support for central gov variable


data$TrustCentGov_new <- case_when(
  data$TrustCentGov == 1 ~ 5,
  data$TrustCentGov  == 2 ~ 4, 
  data$TrustCentGov == 3 ~ 3, 
  data$TrustCentGov == 4 ~ 2,
  data$TrustCentGov == 5 ~ 1,
  TRUE ~ NA_real_
)

# Make support for provincial government variable

data$TrustProvGov_new <- case_when(
  data$TrustProvGov == 1 ~ 5,
  data$TrustProvGov  == 2 ~ 4, 
  data$TrustProvGov == 3 ~ 3, 
  data$TrustProvGov == 4 ~ 2,
  data$TrustProvGov == 5 ~ 1,
  TRUE ~ NA_real_
)


# Conjoint outcome variables 

data$Choice <- data$choice_dummy
data$Rating <- data$rating

# Treatment Variables

data$TreatGroup <- data$TreatInsc

# Set NAs of TreatGroup to 0

data$TreatGroup[is.na(data$TreatGroup)] <- 0

# Make Variable a dummy. Respondents who answered 1,2 or 3
# is response to the end of the prime / manipulation check
# where 1 is true, 2 is false and 3 is Don't know. All of these
# respondents still received the insecurity prime, so we set treatment as any
# answer.

data$TreatGroup <- ifelse(data$TreatGroup > 0, 1, 0)


# Make treatment variable factorial 

data$TreatGroup2 <- ifelse(data$TreatGroup == 1, "Insecurity Prime", "Neutral Text")

# Set levels 

data$TreatGroup2 <- factor(data$TreatGroup2, 
                    levels = c("Neutral Text", "Insecurity Prime"))

# Make a Variable for Treatment and Gender

data$Respondent_Gender <- ifelse(data$Res_Gender == 1, "Male",
                          "Female")

data$TreatGender <- ifelse(data$TreatGroup == 1 & data$Respondent_Gender == "Female",
                      "Insecurity Prime - Female",
                      ifelse(data$TreatGroup == 1 & data$Respondent_Gender == "Male",
                        "Insecurity Prime - Male",     
                        ifelse(data$TreatGroup == 0 & data$Respondent_Gender == "Male",
                         "Neutral Text - Male", "Neutral Text - Female")))       


data$TreatGender <- factor(data$TreatGender,
                    levels = c("Neutral Text - Female", "Neutral Text - Male",
                     "Insecurity Prime - Female", "Insecurity Prime - Male"))


# Respondent is a dove if they want international forces
# to leave (below the median of 2)

data$dove <- ifelse(data$IntForces_new < 3, "Dove", "Hawk")


data$TreatDove <- ifelse(data$TreatGroup == 1 & data$dove == "Dove",
                          "Insecurity Prime - Dove",
                           ifelse(data$TreatGroup == 1 & data$dove == "Hawk",
                                  "Insecurity Prime - Hawk",     
                                  ifelse(data$TreatGroup == 0 & data$dove == "Dove",
                                         "Neutral Text - Dove", "Neutral Text - Hawk")))       


data$TreatDove <- factor(data$TreatDove,
                    levels = c("Neutral Text - Dove", "Neutral Text - Hawk",
                     "Insecurity Prime - Dove", "Insecurity Prime - Hawk"))


# Respondent has high levels of trust in their provincial gov if their trust value is greater
# than three

data$high_prov_trust <- ifelse(data$TrustProvGov_new > 3, 1, 0)

data$high_prov_trust_factor <- ifelse(data$high_prov_trust == 1, "High Trust",
                                      "Low Trust")
data$high_prov_trust_factor <- as.factor(data$high_prov_trust_factor)


# Set baselines for MM analysis
baselines <- list()
baselines$Gender <- "Male"


# Load cregg package fro MM analysis
# install.packages("cregg")
library("cregg") 
# note: cregg masks amce function from cjoint, 
# detach this package before attempting to 
# re-run amce analysis above or call cjoint::amce


######## 3. H1: Insecurity and Support fro Female Leadership #################
##    Calculate marginal means, do F-tests, make
##    tables and plots 


# Calculate marginal means. h0 = 0.5 because it is a forced choice
# TreatGroup2 is Neutral Text vs. Insecurity Prime

stacked_h1 <- cj(data, Choice  ~ Gender,
                 id = ~ ID,
                 estimate = "mm",
                 by = ~ TreatGroup2,
                 h0 = 0.5
)
 


# F-Test results reported in text with choice
f.test.h1 <- cj_anova(data, Choice  ~ Gender,
                       by = ~ TreatGroup2
)


f.test.h1 # p = 0.70 (F-test Result Section 4.1)


# F-Test results reported in text with rating variable
f.test.h1_rating <- cj_anova(data, Rating  ~ Gender,
                      by = ~ TreatGroup2
)


f.test.h1_rating # p = 0.28 (F-test Result Section 4.1)


stacked_h1_rating <- cj(data, Rating  ~ Gender,
                 id = ~ ID,
                 estimate = "mm",
                 by = ~ TreatGroup2,
                 h0 = 2.5
)


# SI Table 7 (SI Section 2.4)

h1 <- stacked_h1[, c("BY", "level", "estimate", "std.error")]

estimate <- h1[, "estimate"]

ses <- h1[, "std.error"]

table <- rbind(estimate, ses)

h1_female <- table[, c(1,3)]
h1_male <- table[,c(2,4)]

h1_all <- rbind(h1_female, h1_male)

colnames(h1_all) <- c("Neutral Text", "Insecurity Prime")

h1_all <- as.table(h1_all)

h1_all <- round(h1_all, digits = 3) 

Gender <- c("Female", "", "Male", "")

h1_choice <- as.data.frame(cbind(Gender, h1_all))

colnames(h1_choice) <- c("Leader's Gender", "Neutral Text", "Insecurity Prime")

view(h1_choice)

# calculation number of observations with 
table(data$TreatGroup2)

# Hand code Table 7 (SI 2.4)


## SI Table 8 H1: Rating (SI 2.4)

h1 <- stacked_h1_rating[, c("BY", "level", "estimate", "std.error")]

estimate <- h1[, "estimate"]

ses <- h1[, "std.error"]

table <- rbind(estimate, ses)

h1_female <- table[, c(1,3)]
h1_male <- table[,c(2,4)]

h1_all <- rbind(h1_female, h1_male)

colnames(h1_all) <- c("Neutral Text", "Insecurity Prime")

h1_all <- as.table(h1_all)

h1_all <- round(h1_all, digits = 3) 

Gender <- c("Female", "", "Male", "")

h1_rating <- as.data.frame(cbind(Gender, h1_all))

colnames(h1_rating) <- c("Leader's Gender", "Neutral Text", "Insecurity Prime")

view(h1_rating)

# calculation number of observations with 
table(data$TreatGroup2)


# Hand code Table 8 into SI


# Figure 1 (Section 4.1 in manuscript)

# Rename TreatGroup2 to Treatment to facilitate reading the 
# graph

data$Treatment <- data$TreatGroup2


h1 <- cj(data, Choice ~ Gender, id = ~ID, estimate = "mm", by = ~ Treatment)

# Make treatment conditions more intuitive 

group <- c("Female Leader Neutral Text", "Male Leader Neutral Text", "Female Leader Insecurity Prime",
           "Male Leader Insecurity Prime")

h1 <- cbind(h1, group)

h1$group2 <- factor(h1$group, levels = 
                     c("Female Leader Neutral Text", "Male Leader Neutral Text", "Female Leader Insecurity Prime",
                     "Male Leader Insecurity Prime"))
                     
            
figure1 <- ggplot(data = as.data.frame(h1),
          aes(x = estimate,
           y = group2)) + 
          geom_point(size = 3)+ 
          geom_pointrange(aes(xmin= estimate - 2*std.error,
           xmax=estimate + 2*std.error))+
          labs(x="Marginal Mean", y = "")+
          geom_vline(xintercept = 0.5, linetype = "dashed")+
          scale_y_discrete(labels=c("Male Leader Neutral Text" = 
          "Male Leader\n(Neutral Text)",
          "Female Leader Neutral Text" = 
            "Female Leader\n(Neutral Text)",
          "Male Leader Insecurity Prime" = 
          "Male Leader\n(Insecurity Prime)",
          "Female Leader Insecurity Prime" = 
          "Female Leader\n(Insecurity Prime)"))

# export Figure 1


## H1: Power Analysis (SI Section 2.13)

# Power Analysis - calculate the Average Marginal Conditional Effect (AMCE) for
# treatment and male leadership attribute.  

amces_choice <- cj(data, Choice ~ Gender, 
                   id = ~ID, estimate = "amce", 
                   by = ~TreatGroup2)


# Insecurity prime has an amce estimate of 0.07 for profiles with Males 

# Use cjpowr package
# to calculate power for interactive effect of 
# insecurity treatment and male leadership attribute
# delta3 (amces) is 0.07. N is 7794. Treat probability
# is 0.25. Four treatment conditions, each likely to be had.

cjpowr_amcie(
  delta1 = 0,
  delta3 = 0.07,
  alpha = 0.05,
  n = 7794,
  levels1 = 2,
  levels2 = 2,
  p00 = 0.25,
  p10= 0.25,
  p01=0.25,
  p11=0.25,
  sims = 1e+05,
  delta0 = 0.5
)

# power estimate of 0.875 (SI Section 2.3)

cjpowr_amcie(
  delta1 = 0,
  delta3 = 0.07,
  alpha = 0.05,
  power = 0.8,
  levels1 = 2,
  levels2 = 2,
  p00 = 0.25,
  p10= 0.25,
  p01=0.25,
  p11=0.25,
  sims = 1e+05,
  delta0 = 0.5
)

# For a power of 0.8, need 6328 observations
# below our n of 7794 


# SI Figure 6: H1 Power Analysis (SI 2.13)

cjpowr_amcie_vec <- Vectorize(cjpowr_amcie)

d <- expand.grid(delta3 = 0.07, n = seq(from = 100, to = 10000, length.out = 1000))

df <- t(cjpowr_amcie_vec(delta3 = d$delta3, n = d$n, sims = 10000, levels1 = 2, levels2=2, 
                         alpha = 0.05, delta0 = 0.5))

df <- data.frame(df)

df$power <- as.numeric(df$power)
df$n <- as.numeric(df$n)


si_figure6 <- ggplot(df, aes(y =power, x=n))+
  geom_line()+ ylim(0,1) +
  geom_vline(xintercept = 7794, linetype = 3, 
             size = 1, color = "black") +
  geom_hline(yintercept = 0.8, linetype = 3, 
             size = 1, color = "black") +
  labs(x = "Sample Size",
       y = "Power")


## H1: Results Comparing Doves and Hawks (SI Section 2.5) ##

stacked_h1b <- cj(data, Choice  ~ Gender,
                 id = ~ ID,
                 estimate = "mm",
                 by = ~ TreatDove,
                 h = 0.5
)

## F-Test results reported in text
f.test.h1b <- cj_anova(data, Choice  ~ Gender,
                      by = ~ TreatDove
)

f.test.h1b # p = 0.99 (Result in SI Section 2.5).


## SI Figure 2 (SI Section 2.5)

# Relevel TreatDove variable

data$TreatDove <- factor(data$TreatDove,
                         levels = c("Neutral Text - Dove",
                                    "Insecurity Prime - Dove",
                                    "Neutral Text - Hawk",
                                    "Insecurity Prime - Hawk"))

h1 <- cj(data, Choice ~ Gender, id = ~ID, estimate = "mm", 
         by = ~ TreatDove)

group <- c("Female Leader Neutral Text - Dove", 
           "Male Leader Neutral Text - Dove",
           "Female Leader Insecurity Prime - Dove",
           "Male Leader Insecurity Prime - Dove",
           "Female Leader Neutral Text - Hawk", 
           "Male Leader Neutral Text - Hawk",
           "Female Leader Insecurity Prime - Hawk",
           "Male Leader Insecurity Prime - Hawk")


h1 <- cbind(h1, group)

# Relevel the treatment  groups 
h1$group2 <- factor(h1$group,
                    levels = c("Female Leader Insecurity Prime - Dove",
                               "Male Leader Insecurity Prime - Dove", 
                               "Female Leader Neutral Text - Dove", 
                               "Male Leader Neutral Text - Dove",
                               "Female Leader Insecurity Prime - Hawk",
                               "Male Leader Insecurity Prime - Hawk", 
                               "Female Leader Neutral Text - Hawk", 
                               "Male Leader Neutral Text - Hawk"
                    ))

# Make figure 2 (in SI)
si_figure2 <- ggplot(data = as.data.frame(h1),
                aes(x = estimate,
                    y = group2)) + 
  geom_point(size = 3)+ 
  geom_pointrange(aes(xmin= estimate - 2*std.error,
                      xmax=estimate + 2*std.error))+
  labs(x="Marginal Mean", y = "")+
  geom_vline(xintercept = 0.5, linetype = "dashed")+
  scale_y_discrete(labels=c("Male Leader Neutral Text - Dove" = 
                              "Male Leader\n(Neutral Text - \n Dove)",
                            "Male Leader Insecurity Prime - Dove" = 
                              "Male Leader\n(Insecurity Prime - \n Dove)",
                            "Female Leader Neutral Text - Dove" = 
                              "Female Leader\n(Neutral Text - \n Dove)",
                            "Female Leader Insecurity Prime - Dove" = 
                              "Female Leader\n(Insecurity Prime - \n Dove)",
                            "Male Leader Neutral Text - Hawk" = 
                              "Male Leader\n(Neutral Text - \n Hawk)",
                            "Male Leader Insecurity Prime - Hawk" = 
                              "Male Leader\n(Insecurity Prime - \n Hawk)",
                            "Female Leader Neutral Text - Hawk" = 
                              "Female Leader\n(Neutral Text - \n Hawk)",
                            "Female Leader Insecurity Prime - Hawk" = 
                              "Female Leader\n(Insecurity Prime - \n Hawk)"))

# Export SI Figure 2 (SI Section 2.5)


##### 4. H2: Insecurity and Support for Female Leadership Across Genders #################################
##    Test insecurity prime's effect
##    across genders. Calculate marginal means,
##    do F tests Create plots and tables
##    for male and female subsets. 
##    Run power analyses


## H2 Insecurity prime will decrease men's preferences for 
## female leaders


# Subset to male respondents 
male <- data[data$Respondent_Gender == "Male",]


# Calculate marginal means
stacked_h2 <- cj(male, Choice  ~ Gender,
                  id = ~ ID,
                  estimate = "mm",
                  by = ~ TreatGroup2,
                  h = 0.5
)


# F-Test results reported in text
f.test.h2 <- cj_anova(male, Choice  ~ Gender,
                       by = ~ TreatGroup2
) 

f.test.h2 # f test is 0.485 (in Manuscript Section 4.2)
# Reject H2b


## SI Table 9: H2: Marginal Means Preferences for Female Leadership Across
## Gender and Treatment Groups (SI Section 2.6)

# Calculate marginal means


stacked_h2 <- cj(data, Choice  ~ Gender,
                  id = ~ ID,
                  estimate = "mm",
                  by = ~ TreatGender,
                  h = 0.5
)


h2 <- stacked_h2[, c("BY", "level", "estimate", "std.error", "TreatGender")]

estimate <- h2[, "estimate"]

ses <- h2[, "std.error"]

table <- rbind(estimate, ses)

## Female Preferences

table_female <- table[, c(1:2,5:6)]

table_female_1 <- table_female[, c(1,3)]
table_female_2 <- table_female[, c(2,4)]

table_female <- rbind(table_female_1, table_female_2)

h2_female <- round(table_female, digits = 3) 

Gender <- c("Female", "", "Male", "")

h2_female <- as.data.frame(cbind(Gender, h2_female))

colnames(h2_female) <- c("Leader's Gender", "Neutral Text", "Insecurity Prime")

view(h2_female) # Handcode into SI

## Male Preferences

table_male <- table[, c(3:4,7:8)]

table_male_1 <- table_male[, c(1,3)]
table_male_2 <- table_male[, c(2,4)]

table_male <- rbind(table_male_1, table_male_2)

h2_male <- round(table_male, digits = 3) 

Gender <- c("Female", "", "Male", "")

h2_male <- as.data.frame(cbind(Gender, h2_male))

colnames(h2_male) <- c("Leader's Gender", "Neutral Text", "Insecurity Prime")

view(h2_male) #handcode into SI


## H2: Figure 2 (Manuscript Section 4.2)

# Relevel TreatGender variable

data$TreatGender <- factor(data$TreatGender,
                  levels = c("Neutral Text - Female",
                             "Insecurity Prime - Female",
                             "Neutral Text - Male",
                             "Insecurity Prime - Male"))
                             
h2 <- cj(data, Choice ~ Gender, id = ~ID, estimate = "mm", 
         by = ~ TreatGender)

group <- c("Female Leader Neutral Text - Female Respondent", 
           "Male Leader Neutral Text - Female Respondent",
           "Female Leader Insecurity Prime - Female Respondent",
           "Male Leader Insecurity Prime - Female Respondent",
           "Female Leader Neutral Text - Male Respondent", 
           "Male Leader Neutral Text - Male Respondent",
           "Female Leader Insecurity Prime - Male Respondent",
           "Male Leader Insecurity Prime - Male Respondent")
           
           
h2 <- cbind(h2, group)

# Relevel the treatment  groups 
h2$group2 <- factor(h2$group,
              levels = c("Female Leader Insecurity Prime - Female Respondent",
                         "Male Leader Insecurity Prime - Female Respondent", 
                         "Female Leader Neutral Text - Female Respondent", 
                              "Male Leader Neutral Text - Female Respondent",
                              "Female Leader Insecurity Prime - Male Respondent",
                         "Male Leader Insecurity Prime - Male Respondent", 
                             "Female Leader Neutral Text - Male Respondent", 
                              "Male Leader Neutral Text - Male Respondent"
                              ))


## Make plot for only male respondents

h2_male_a <- h2[h2$group2 == "Female Leader Neutral Text - Male Respondent",]
h2_male_b <- h2[h2$group2 == "Male Leader Neutral Text - Male Respondent",]
h2_male_c <- h2[h2$group2 == "Female Leader Insecurity Prime - Male Respondent",]
h2_male_d <- h2[h2$group2 == "Male Leader Insecurity Prime - Male Respondent",]

h2_male <- rbind(h2_male_a, h2_male_b, h2_male_c, h2_male_d)


## H2: Figure 2 (Male Respondents Only) (in Manuscript)
                
figure2 <- ggplot(data = as.data.frame(h2_male),
                aes(x = estimate,
                    y = group2)) + 
  geom_point(size = 3)+ 
  geom_pointrange(aes(xmin= estimate - 2*std.error,
                      xmax=estimate + 2*std.error))+
  labs(x="Marginal Mean", y = "")+
  geom_vline(xintercept = 0.5, linetype = "dashed")+
  scale_y_discrete(labels=c("Male Leader Neutral Text - Male Respondent" = 
                              "Male Leader\n(Neutral Text)",
                            "Male Leader Insecurity Prime - Male Respondent" = 
                              "Male Leader\n(Insecurity Prime)",
                            "Female Leader Neutral Text - Male Respondent" = 
                              "Female Leader\n(Neutral Text)",
                            "Female Leader Insecurity Prime - Male Respondent" = 
                              "Female Leader\n(Insecurity Prime)"))

# Export Figure 2

## H2: Hawks vs. Doves (Men only)

stacked_h2_male_doves <- cj(male, Choice  ~ Gender,
                 id = ~ ID,
                 estimate = "mm",
                 by = ~ TreatDove,
                 h = 0.5
)

## F-Test results reported in text
f.test.h2_males_doves <- cj_anova(male, Choice  ~ Gender,
                      by = ~ TreatDove
)

f.test.h2_males_doves # p = 0.73 (Section 4.2 in Manuscript)


## SI Figure 3

h2_male_dove <- cj(male, Choice ~ Gender, id = ~ID, estimate = "mm", 
          by = ~ TreatDove)

group <- c("Female Leader Neutral Text - Dove", 
           "Male Leader Neutral Text - Dove",
           "Female Leader Insecurity Prime - Dove",
           "Male Leader Insecurity Prime - Dove",
           "Female Leader Neutral Text - Hawk", 
           "Male Leader Neutral Text - Hawk",
           "Female Leader Insecurity Prime - Hawk",
           "Male Leader Insecurity Prime - Hawk")


h2_male_dove  <- cbind(h2_male_dove, group)

# Relevel the treatment  groups 
h2_male_dove$group2 <- factor(h2_male_dove$group,
                     levels = c("Female Leader Insecurity Prime - Dove",
                                "Male Leader Insecurity Prime - Dove", 
                                "Female Leader Neutral Text - Dove", 
                                "Male Leader Neutral Text - Dove",
                                "Female Leader Insecurity Prime - Hawk",
                                "Male Leader Insecurity Prime - Hawk", 
                                "Female Leader Neutral Text - Hawk", 
                                "Male Leader Neutral Text - Hawk"
                     ))

# Make SI Figure 3 (SI Section 2.7)

si_figure3 <- ggplot(data = as.data.frame(h2_male_dove),
                  aes(x = estimate,
                      y = group2)) + 
  geom_point(size = 3)+ 
  geom_pointrange(aes(xmin= estimate - 2*std.error,
                      xmax=estimate + 2*std.error))+
  labs(x="Marginal Mean", y = "")+
  geom_vline(xintercept = 0.5, linetype = "dashed")+
  scale_y_discrete(labels=c("Male Leader Neutral Text - Dove" = 
                              "Male Leader\n(Neutral Text - \n Dove)",
                            "Male Leader Insecurity Prime - Dove" = 
                              "Male Leader\n(Insecurity Prime - \n Dove)",
                            "Female Leader Neutral Text - Dove" = 
                              "Female Leader\n(Neutral Text - \n Dove)",
                            "Female Leader Insecurity Prime - Dove" = 
                              "Female Leader\n(Insecurity Prime - \n Dove)",
                            "Male Leader Neutral Text - Hawk" = 
                              "Male Leader\n(Neutral Text - \n Hawk)",
                            "Male Leader Insecurity Prime - Hawk" = 
                              "Male Leader\n(Insecurity Prime - \n Hawk)",
                            "Female Leader Neutral Text - Hawk" = 
                              "Female Leader\n(Neutral Text - \n Hawk)",
                            "Female Leader Insecurity Prime - Hawk" = 
                              "Female Leader\n(Insecurity Prime - \n Hawk)"))


## export  SI figure 3


## H2: Insecurity and Support for Female Leadership Among Female Respondents

## H2: Female Respondents (Manuscript Section 4.2)

# Subset data to only female respondents

female <- data[data$Respondent_Gender == "Female",]


# Calculate marginal means
stacked_h2 <- cj(female, Choice  ~ Gender,
                 id = ~ ID,
                 estimate = "mm",
                 by = ~ TreatGroup2,
                 h = 0.5
)


# F-Test results reported in text
f.test.h2 <- cj_anova(female, Choice  ~ Gender,
                      by = ~ TreatGroup2
) 

f.test.h2 # F test is 0.046 (Results in Manuscript Section 4.2)



# F-Test results reported in text
f.test.h2c_rating <- cj_anova(female, Rating  ~ Gender,
                              by = ~ TreatGroup2
) 

f.test.h2c_rating# f test is 0.04



## Figure 3 in Manuscript

h2_female_a <- h2[h2$group2 == "Female Leader Neutral Text - Female Respondent",]
h2_female_b <- h2[h2$group2 == "Male Leader Neutral Text - Female Respondent",]
h2_female_c <- h2[h2$group2 == "Female Leader Insecurity Prime - Female Respondent",]
h2_female_d <- h2[h2$group2 == "Male Leader Insecurity Prime - Female Respondent",]

h2_female <- rbind(h2_female_a, h2_female_b, h2_female_c, h2_female_d)


figure3 <- ggplot(data = as.data.frame(h2_female),
                aes(x = estimate,
                    y = group2)) + 
  geom_point(size = 3)+ 
  geom_pointrange(aes(xmin= estimate - 2*std.error,
                      xmax=estimate + 2*std.error))+
  labs(x="Marginal Mean", y = "")+
  geom_vline(xintercept = 0.5, linetype = "dashed")+
  scale_y_discrete(labels=c("Male Leader Neutral Text - Female Respondent" = 
                              "Male Leader\n(Neutral Text)",
                            "Male Leader Insecurity Prime - Female Respondent" = 
                              "Male Leader\n(Insecurity Prime)",
                            "Female Leader Neutral Text - Female Respondent" = 
                              "Female Leader\n(Neutral Text)",
                            "Female Leader Insecurity Prime - Female Respondent" = 
                              "Female Leader\n(Insecurity Prime)"))



## Examine Female Hawks and Doves


stacked_h2_female_doves <- cj(female, Choice  ~ Gender,
                 id = ~ ID,
                 estimate = "mm",
                 by = ~ TreatDove,
                 h = 0.5
)

## F-Test results reported in text
f.test.h2_female_doves <- cj_anova(female, Choice  ~ Gender,
                      by = ~ TreatDove
)

f.test.h2_female_doves # p = 0.14 (Reported in manuscript).


# SI Figure 4


h2_female_doves <- cj(female, Choice ~ Gender, id = ~ID, estimate = "mm", 
          by = ~ TreatDove)

group <- c("Female Leader Neutral Text - Dove", 
           "Male Leader Neutral Text - Dove",
           "Female Leader Insecurity Prime - Dove",
           "Male Leader Insecurity Prime - Dove",
           "Female Leader Neutral Text - Hawk", 
           "Male Leader Neutral Text - Hawk",
           "Female Leader Insecurity Prime - Hawk",
           "Male Leader Insecurity Prime - Hawk")


h2_female_doves <- cbind(h2_female_doves, group)

# Relevel the treatment  groups 
h2_female_doves$group2 <- factor(h2_female_doves$group,
                     levels = c("Female Leader Insecurity Prime - Dove",
                                "Male Leader Insecurity Prime - Dove", 
                                "Female Leader Neutral Text - Dove", 
                                "Male Leader Neutral Text - Dove",
                                "Female Leader Insecurity Prime - Hawk",
                                "Male Leader Insecurity Prime - Hawk", 
                                "Female Leader Neutral Text - Hawk", 
                                "Male Leader Neutral Text - Hawk"
                     ))

# SI Figure 4 (SI Section 2.7)

si_figure4 <- ggplot(data = as.data.frame(h2_female_doves),
                   aes(x = estimate,
                       y = group2)) + 
  geom_point(size = 3)+ 
  geom_pointrange(aes(xmin= estimate - 2*std.error,
                      xmax=estimate + 2*std.error))+
  labs(x="Marginal Mean", y = "")+
  geom_vline(xintercept = 0.5, linetype = "dashed")+
  scale_y_discrete(labels=c("Male Leader Neutral Text - Dove" = 
                              "Male Leader\n(Neutral Text - \n Dove)",
                            "Male Leader Insecurity Prime - Dove" = 
                              "Male Leader\n(Insecurity Prime - \n Dove)",
                            "Female Leader Neutral Text - Dove" = 
                              "Female Leader\n(Neutral Text - \n Dove)",
                            "Female Leader Insecurity Prime - Dove" = 
                              "Female Leader\n(Insecurity Prime - \n Dove)",
                            "Male Leader Neutral Text - Hawk" = 
                              "Male Leader\n(Neutral Text - \n Hawk)",
                            "Male Leader Insecurity Prime - Hawk" = 
                              "Male Leader\n(Insecurity Prime - \n Hawk)",
                            "Female Leader Neutral Text - Hawk" = 
                              "Female Leader\n(Neutral Text - \n Hawk)",
                            "Female Leader Insecurity Prime - Hawk" = 
                              "Female Leader\n(Insecurity Prime - \n Hawk)"))


# export to SI


#### 5. H2: Robustness Checks #########


# SI Section 2.9 Exclude highly educated (any university degree)

edu <- data[data$edu_numeric < 5,]

h2_edu <- cj(edu, Choice ~ Gender, id = ~ID, estimate = "mm", 
             by = ~ TreatGender)

# F-Test results reported in text
f.test.h2_edu <- cj_anova(edu, Choice  ~ Gender,
                          by = ~ TreatGender
)

f.test.h2_edu # stat difference at 0.001


# SI Table 11 

edu_table <- h2_edu[, c("BY", "level", "estimate", "std.error", "TreatGender")]

estimate <- edu_table[, "estimate"]

ses <- edu_table[, "std.error"]

table <- rbind(estimate, ses)

## Female Preferences


table_female <- table[, c(1:4)]

table_female_1 <- table_female[, c(1,3)]
table_female_2 <- table_female[, c(2,4)]

table_female <- rbind(table_female_1, table_female_2)

h2_female <- round(table_female, digits = 3) 

Gender <- c("Female", "", "Male", "")

h2_female <- as.data.frame(cbind(Gender, h2_female))

colnames(h2_female) <- c("Leader's Gender",  "Neutral Text", "Insecurity Prime")

view(h2_female)

## Male Table

table_male <- table[, c(5:8)]

table_male_1 <- table_male[, c(1,3)]
table_male_2 <- table_male[, c(2,4)]

table_male <- rbind(table_male_1, table_male_2)

h2_male <- round(table_male, digits = 3) 

Gender <- c("Female", "", "Male", "")

h2_male <- as.data.frame(cbind(Gender, h2_male))

colnames(h2_male) <- c("Leader's Gender", "Neutral Text", "Insecurity Prime")

view(h2_male)

# find observations

table(edu$TreatGender)


# Handcode Table 11 in the SI


edu_female <- edu[edu$Respondent_Gender == "Female",]

h2_edu_female <- cj(edu_female, Choice ~ Gender, id = ~ID, estimate = "mm", 
                    by = ~ TreatGroup2)


## F-Test results reported in text
f.test.edu_female <- cj_anova(edu_female, Choice  ~ Gender,
                              by = ~ TreatGroup2
)


f.test.edu_female # stat difference
# at 0.04


# Less educated men

edu_male <- edu[edu$Respondent_Gender == "Male",]

h2_edu_male <- cj(edu_male, Choice ~ Gender, id = ~ID, estimate = "mm", 
                  by = ~ TreatGroup2)


## F-Test results reported in text
f.test.edu_male <- cj_anova(edu_male, Choice  ~ Gender,
                            by = ~ TreatGroup2
)


f.test.edu_male # stat difference
# at 0.49


# SI Section 2.10 Respondents outside of Kunduz


no_kundoz <- data[data$Province != "Kundoz",]

h2_no_kundoz <- cj(no_kundoz, Choice ~ Gender, id = ~ID, estimate = "mm", 
                   by = ~ TreatGender)



## SI Table 12

h2_no_kundoz <- h2_no_kundoz[, c("BY", "level", "estimate", "std.error", "TreatGender")]

estimate <- h2_no_kundoz[, "estimate"]

ses <- h2_no_kundoz[, "std.error"]

table <- rbind(estimate, ses)

# Female Preferences

table_female <- table[, c(1:4)]

table_female_1 <- table_female[, c(1,3)]
table_female_2 <- table_female[, c(2,4)]

table_female <- rbind(table_female_1, table_female_2)

h2_female <- round(table_female, digits = 3) 

Gender <- c("Female", "", "Male", "")

h2_female <- as.data.frame(cbind(Gender, h2_female))

colnames(h2_female) <- c("Leader's Gender",  "Neutral Text", "Insecurity Prime")

view(h2_female)


# Male Preferences

table_male <- table[, c(5:8)]

table_male_1 <- table_male[, c(1,3)]
table_male_2 <- table_male[, c(2,4)]

table_male <- rbind(table_male_1, table_male_2)

h2_male <- round(table_male, digits = 3) 

Gender <- c("Female", "", "Male", "")

h2_male <- as.data.frame(cbind(Gender, h2_male))

colnames(h2_male) <- c("Leader's Gender", "Neutral Text", "Insecurity Prime")

view(h2_male)

# find observations

table(no_kundoz$TreatGender)

## Handcode table in SI

# h2: subset to women who don't live in Kunduz

no_kundoz_female <- no_kundoz[no_kundoz$Respondent_Gender == "Female",]

h2_no_kundoz_female <- cj(no_kundoz_female, Choice ~ Gender, id = ~ID, estimate = "mm", 
                          by = ~ TreatGroup2)


## F-Test results reported in SI 2.10
f.test.h2_no_kundoz_no_female <- cj_anova(no_kundoz_female, Choice  ~ Gender,
                                          by = ~ TreatGroup2
) 

f.test.h2_no_kundoz_no_female # f test p: 0.01


# Subset to males who don't live in kundoz

no_kundoz_male <- no_kundoz[no_kundoz$Respondent_Gender == "Male",]

h2_no_kundoz_male <- cj(no_kundoz_male, Choice ~ Gender, id = ~ID, estimate = "mm", 
                        by = ~ TreatGroup2)


# F-Test results reported in text
f.test.h2_no_kundoz_male <- cj_anova(no_kundoz_male, Choice  ~ Gender,
                                     by = ~ TreatGroup2
) 


# F-test reported in SI 2.10
f.test.h2_no_kundoz_male # f.test: p: 0.968


# SI Section 2.11 Unemployed Respondents


emp <- data[data$Employ_Unem == 1,]


h2_emp <- cj(emp, Choice ~ Gender, id = ~ID, estimate = "mm", 
             by = ~ TreatGender)


## SI Table 13 

emp_table <- h2_emp[, c("BY", "level", "estimate", "std.error", "TreatGender")]

estimate <- emp_table[, "estimate"]

ses <- emp_table[, "std.error"]

table <- rbind(estimate, ses)

## Female Preferences

table_female <- table[, c(1:4)]

table_female_1 <- table_female[, c(1,3)]
table_female_2 <- table_female[, c(2,4)]

table_female <- rbind(table_female_1, table_female_2)

h2_female <- round(table_female, digits = 3) 

Gender <- c("Female", "", "Male", "")

h2_female <- as.data.frame(cbind(Gender, h2_female))

colnames(h2_female) <- c("Leader's Gender",  "Neutral Text", "Insecurity Prime")

view(h2_female)


## Male Preferences

table_male <- table[, c(5:8)]

table_male_1 <- table_male[, c(1,3)]
table_male_2 <- table_male[, c(2,4)]

table_male <- rbind(table_male_1, table_male_2)

h2_male <- round(table_male, digits = 3) 

Gender <- c("Female", "", "Male", "")

h2_male <- as.data.frame(cbind(Gender, h2_male))

colnames(h2_male) <- c("Leader's Gender", "Neutral Text", "Insecurity Prime")

view(h2_male)

# find observations

table(emp$TreatGender)

# Hand code table into text


# Subset to Unemployed Women


emp_female <- emp[emp$Respondent_Gender == "Female",]


h2_emp_female <- cj(emp_female, Choice ~ Gender, id = ~ID, estimate = "mm", 
                    by = ~ TreatGroup2)


# F-Test results reported in text
f.test.emp_female <- cj_anova(emp_female, Choice  ~ Gender,
                              by = ~ TreatGroup2
)


# Reported in SI Section 2.11
f.test.emp_female # f test: 0.035

# Male unemployed

emp_male <- emp[emp$Respondent_Gender == "Male",]


h2_emp_male <- cj(emp_male, Choice ~ Gender, id = ~ID, estimate = "mm", 
                  by = ~ TreatGroup2)


# F-Test results reported in text
f.test.emp_male <- cj_anova(emp_male, Choice  ~ Gender,
                            by = ~ TreatGroup2
)


# Reported in SI Section 2.11
f.test.emp_male #p.value = 0.93


## H2 Power Analysis

## Power Analysis for H2

# SI Figure 7: H2 Power Analysis (Male Respondents Only) (SI Section 2.13)

amces_choice_h2 <- cj(male, Choice ~ Gender, 
                      id = ~ID, estimate = "amce", 
                      by = ~TreatGroup2)


# Estimate an AMCIE for treatment and male leadership
# attribute of 0.08. There are 3,636 male respondent 
# profile observations 

cjpowr_amcie(
  delta1 = 0,
  delta3 = 0.08,
  alpha = 0.05,
  n = 3636,
  levels1 = 2,
  levels2 = 2,
  p00 = 0.25,
  p01=0.25,
  p10=0.25,
  p11=0.25,
  sims = 1e+05,
  delta0 = 0.5
)


# power: 0.68

cjpowr_amcie(
  delta1 = 0,
  delta3 = 0.08,
  alpha = 0.05,
  power = 0.8,
  levels1 = 2,
  levels2 = 2,
  p00 = 0.25,
  p01=0.25,
  p10=0.25,
  p11=0.25,
  sims = 1e+05,
  delta0 = 0.5
)

# need 4827 profile observations to have power of 0.8. 
# this is more than current number of male observations: 3636

# make plot 

cjpowr_amcie_vec <- Vectorize(cjpowr_amcie)


d <- expand.grid(delta3 = 0.08, n = seq(from = 100, to = 10000, length.out = 1000))

df <- t(cjpowr_amcie_vec(delta3 = d$delta3, n = d$n, sims = 10000, levels1 = 2, levels2=2, 
                         alpha = 0.05, delta0 = 0.5))

df <- data.frame(df)

df$power <- as.numeric(df$power)
df$n <- as.numeric(df$n)

## Make plot with power as the y axis, and n the x axis

si_figure7 <- ggplot(df, aes(y =power, x=n))+
  geom_line()+ ylim(0,1) +
  geom_vline(xintercept = 3636, linetype = 3, 
             size = 1, color = "black") +
  geom_hline(yintercept = 0.8, linetype = 3, 
             size = 1, color = "black") +
  labs(x = "Sample Size",
       y = "Power")


## H2 power analysis for female respondents

# H2: Power Analysis (Female Respondents)


amces_choice_h2 <- cj(female, Choice ~ Gender, 
                      id = ~ID, estimate = "amce", 
                      by = ~TreatGroup2)


# Estimate an AMCIE for treatment and male leadership
# attribute of 0.06. There are 4,158 female respondents 
# profile observations 

cjpowr_amcie(
  delta1 = 0,
  delta3 = 0.06,
  alpha = 0.05,
  n = 4158,
  levels1 = 2,
  levels2 = 2,
  p00 = 0.25,
  p01=0.25,
  p10=0.25,
  p11=0.25,
  sims = 1e+05,
  delta0 = 0.5
)


# power: 0.493, less than 0.8 

cjpowr_amcie(
  delta1 = 0,
  delta3 = 0.06,
  alpha = 0.05,
  power = 0.8,
  levels1 = 2,
  levels2 = 2,
  p00 = 0.25,
  p01=0.25,
  p10=0.25,
  p11=0.25,
  sims = 1e+05,
  delta0 = 0.5
)

# need 8642 profile observations to have power of 0.8. 
# type_s errors: <0.00, exp_typeM: 1.13








